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Abstract. The time evolution of initially balanced, rapidly rotating models for 
an isolated disk of highly flattened galaxies of stars is calculated. The method of 
direct integration of the Newtonian equations of motion of stars over a time span of 
many galactic rotations is applied. Use of modern concurrent supercomputers has 
enabled us to make long simulation runs using a sufficiently large number of particles 
N — 30, 000. One of the goals of the present simulation is to test the validities of a 
modified local criterion for stability of Jeans-type gravity perturbations (e.g. those 
produced by a barlike structure, a spontaneous perturbation and/or a companion 
galaxy) in a self-gravitating, mfinitesimally thin and collisionless disk. In addition to 
the local criterion we are interested in how model particles diffuse in velocity. This 
is of considerable interest in the kinetic theory of stellar disks. Certain astronomical 
implications of the simulations to actual disk-shaped (i.e. rapidly rotating) galax- 
ies are explored. The weakly nonlinear, or quasi-linear kinetic theory of the Jeans 
instability in disk galaxies of stars is described as well. 
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1. Introduction 

In galactic astronomy, the most intensely studied of the problems of 
pattern formation is the problem of the formation of spiral struc- 
ture, which has brought forth numerous theoretical and numerical ap- 
proaches. Even though at the present time there exists, as yet, no 
satisfactory theory of the origin and conservation of the spiral structure 
so prominent in the Milky Way and many other giant flattened galaxies, 
the majority of experts in the field has yielded to the opinion that the 
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study of the stability of the small-amplitude cooperative oscillations in 
disk-shaped galaxies of stars is the first step towards an understanding 
of the phenomena. This is because in the Milky Way and many other 
disk galaxies, the bulk of the luminous mass, probably > 90%, is in 
the stars. Therefore, it seems likely that the spiral structures are in- 
timately connected with the stellar constituent of a galaxy, and that 
stellar dynamical phenomena play a basic role. Because of the long- 
range nature of the gravitational force between stars, a galaxy exhibits 
collective modes of motions — modes in which the stars in large regions 
move coherently or in unison. A galaxy is characterized also by a high 
"temperature" (a dispersion of random velocities of stars) and hence 
a high "thermal" energy; this thermal energy is much larger than the 
interaction potential between pairs of stars. Because of this, binary 
encounters produce only small-scale scattering of the motions of the 
stars. Over long enough times, these gradually build up to produce large 
deflections that constitute collisions. In many cases these collisional 
effects are so weak that we consider the galaxy as collisionless on the 
Hubble time T ~ 10 10 yr (Chandrasekhar, 1960; Binney and Tremaine, 
1987). Thus, the galaxy is dominated by the collective motions and the 
free streaming of the stars (kinetic effects). 

As usual it is very difficult determine a priori whether a particular 
linearization of nonlinear equations made in analytical studies consti- 
tutes a valid approximation. This can be determined, however, by direct 
computer simulations of the nonlinear theory that mimic the behaviour 
of stars in galaxies. Such modeling gives more detailed information than 
can be obtained analytically, so that the important physical phenomena 
can be determined. In this paper, the results of the approximate theo- 
retical analysis described in the Appendix are compared with our own 
many-body (iV-body) simulations. One of the goals of the simulation is 
to test the validities of a modified local criterion for stability of Jeans- 
type gravity perturbations in a self-gravitating, infinitesimally thin and 
collisionless disk. The fact that the nonaxisymmetric perturbations in 
the differentially rotating system are more unstable than the axisym- 
metric ones is taken into account in this modified Toomre-like (Toomre, 
1964, 1977) criterion. 

In addition to the stability criterion we are interested in how model 
particles diffuse in velocity. As for the present study, such random veloc- 
ity diffusion is caused by stellar disk turbulence. In particular, we verify 
the analytical prediction that as the Jeans-unstable perturbations grow 
the mean-square velocity increases linearly with time. The observation 
of this behavior is a sensitive test of the model as well as the kinetic 
theory of stellar disks. 
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In fact, Morozov (1981) and Griv et al. (1999a) already attempted 
to confirm the modified criterion numerically. However, because of the 
very small number of model stars, ./V = 200, Morozov's results are 
subject to considerable uncertainties, and additional simulations are 
clearly required to settle the issue. Furthermore, for that number of 
model particles, the two-body relaxation time scale is comparable to 
the crossing time, even with Morozov's modest softening parameter, 
raising some question about the applicability of his simulations to 
actual almost collisionless galaxies. Increasing the number of model 
stars is definitely a more reliable procedure. In turn, Griv et al. used a 
different numerical approach of the so-called local iV-body simulations. 
The TV-body experiments in a local or Hill's approximation has been 
pioneered by Toomre (1990), Toomre and Kalnajs (1991) and Salo 
(1995). In these simulations dynamics of particles in small regions of 
the disk are assumed to be statistically independent of dynamics of 
particles in other regions. The local numerical model thus simulates 
only a small part of the system and more distant parts are represented 
as copies of the simulated region. Unlike Morozov (1981) and Griv et al. 
(1999a), we study some aspects of dynamical behavior of stellar systems 
by global simulations using a sufficiently large number of particles. 



2. iV-body simulations 



Different methods are currently employed to simulate global evolution 
of collisionless point-mass systems of galaxies by iV-body experiments. 
One algorithm for a simulation code that has found wide use was de- 
veloped by Hohl and Hockney (1969), Hohl (1971, 1972), Miller (1971), 
Hockney and Brownrigg (1974), Quirk (1972), Athanassoula and Sell- 
wood (1986), Combes et al. (1990), Pfenniger and Friedli (1993), Mer- 
ritt and Sellwood (1994), Donner and Thomasson (1994), and others. 
The method samples the density field with a usually uniform grid. The 
Poisson equation is then solved on this grid using one of a number of 
the fast solvers that are available, usually Fast Fourier Transforms. This 
is an analog of plasma particle-in-cell (PIC) codes. It is believed that 
simulating many billions of stars in actual galaxies by using only sev- 
eral ten or hundred thousands representative stars in PIC experiments 
will be enough to capture the essential physics, which includes the 
ordinary star-star binary relaxation and wave-like collective motion. 
In other fields, such as the simulation of spiral structures, PIC codes 
may be used with moderate success. This is because these fine-scale 
< 1 kpc structures can basically be governed by collisionless processes. 
By increasing the number of cells to reproduce the microstructures, one 
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reduces the average number of particles per cell, and thus increases the 
undesirable effect of particle gravitational encounters. In addition, if 
one uses PIC codes, then one cannot resolve density fluctuations with 
wavelengths smaller than the size of a cell d. This puts a lower limit on 
the wavelengths that it makes any sense to keep; this lower limit is given 
by A m j n ~ d. Other disadvantages of PIC codes are possibilities of late- 
time loss of energy conservation, finite-grid and aliasing instabilities, 
and other numerical artifacts. 

PIC method works well only for regions where there are many 
particles per cell. The problem is serious because there are regions in 
the phase space which are important to the problem but in which the 
distribution function is small. In general, since we will at most have a 
hundred particles in one cell we will have fluctuation of the order of 
10% or more which for many problems is unacceptable. 

In another type of iV-body simulation, Ostriker and Peebles (1973), 
Morozov (1981), Grivnev (1985), Griv et al. (1994), Griv and Zhyt- 
nikov (1995), Griv and Chiueh (1998) investigated a stellar disk by 
direct integration of Newtonian equations of motion. Following them, 
we simulate here the evolution of a model for an isolated thin disk of 
a galaxy by direct integration of the equations of motion of stars: the 
force on a star i due to all other stars is given by 



where G is the constant of gravitation, m 8 , rrij, and rj are the mass, 
position and velocity of the i-th and of the j-th particle, respectively, 
and r cu t is the so-called cutoff radius. This "softening" of the gravita- 
tional potential is a device often used in iV-body simulations to avoid 
numerical difficulties at very close but rare encounters (Miller, 1971; 
Romeo, 1998). The constant parameter r cut reduces the interaction at 
short ranges and puts a lower limit on the size of the particles, i.e. 
the stars in the system can no longer be considered as point-masses - 
they are in fact Plummer spheres with a scale size r cut . In addition, the 
softening parameter r cut leads to a more "collisionless" situation. 

As is known, the latter iV-body models suffer from graininess due 
to the finite number of particles. The graininess manifests itself as a 
numerical "noise" in the calculation. The use of many model particles 
N > 10 4 decreases drastically the noise level. 

Although the direct integration method sounds simple and straigh- 
forward, practical computational limitations require the use of only 
few ten thousands representative stars. This is because we evaluate the 
sum on the right-hand side of Eq. (1) for each star or N times for N 




(1) 
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stars. The sum itself contains N terms; each term requires a number of 
arithmetic operations — for the sake of estimating, let us say, ten. The 
total number of arithmetic operations required to evaluate the force will 
be of the order of 10 x N 2 . For the standard Runge-Kutta calculation 
involving 30, 000 stars the total number of operations would be about 
4 x 10 10 . Assuming that we can do an operation in 10 -8 sec, simply eval- 
uating the forces would require about 6 — 7 min. A typical calculation 
requires several thousand time steps so that ~ 1,000 h would be re- 
quired. Calculations of this magnitude are totally impractical for using 
such models to explore the physics of galaxies. The problem is solved 
by using the nowadays concurrent supercomputers: with the multiple 
processors in supercomputers it is possible to reduce the running time 
by processing more than one group of particles at a time. The use of, 
say, 100 parallel computers would make such calculations practical. In 
the present study, we explore the Inter University 112 processors SGI 
Origin 2000 supercomputer in Israel. 

Note that PIC codes only formally utilized a larger number of par- 
ticles than direct integration codes. This is because in PIC experiments, 
the disk plane was covered with a difference grid. The stars within grid 
cells as well as within a whole model did not interact gravitationally; 
the computation of the mass density in cells followed a fast Fourier 
transform to obtain the potential only at a restricted set of points - 
at the cell centers — by interpolation rules and thus to find the very 
approximate forces at the positions of the particles. This avoids the 
need to compute forces between particle pairs, and makes the amount of 
computation necessary to obtain the forces independent of N. Clearly, 
nothing similar exists in nature, and therefore it follows that one cannot 
compare the number of particles used in these different kinds of sim- 
ulations in any way at all. Modern parallel computing technology on 
concurrent processors had led to a dramatic increase in the computing 
power. The latter allows us to use direct integration methods with a 
sufficiently large number of model stars N = 30, 000. 



3. Numerical model 

In general, the numerical procedure is first to seek stationary solutions 
to the Boltzmann kinetic equation in the self-consistent field approx- 
imation for the density distribution, the angular velocity of regular 
rotation, the dispersion of random velocities of the stars, and then 
to calculate the stability of those solutions to small-amplitude gravity 
perturbations. 
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At the start of the A-body integration, our simulation initializes 
the particles on a set of 100 circular rings with a circular velocity V ro t 
of galactic rotation in the r, ip— plane; the system is isolated in vacuum. 
Consider an uniformly rotating model disk of stars with a surface mass 
density variation given by 



a (r) = a(0)Vl-r 2 AR 2 , (2) 

where c(0) is the central surface density and R is the radius of the 
initial disk. As a solution of a time-independent collisionless Boltzmann 
equation, to ensure initial equilibrium the uniform angular velocity to 
balance the zero-velocity-dispersion disk, 

n = Try/(Mp)/2R, (3) 

was adopted (Hohl, 1971, 1972). Then the position of each particle 
was slightly perturbed by applying a pseudorandom number generator. 
For the uniformly rotating disk, the Maxwellian-distributed random 
velocities with radial c r and azimuthal Cu, dispersions in the plane z = 
according to the familiar Toomre's criterion (Toomre, 1964, 1977) , 



c T 



= = 0.341^07^^, (4) 



may be added (c r = c v ) to the initial circular velocities V TOt = rf^o, 
where k = 20 [1 + (r /VL)(dVL/dr)} 1 / 2 is the ordinary epicyclic frequency 
(Hohl, 1971, 1972). It is crucial to realize that in this case, according 
to Lau and Bertin (1978), Lin and Lau (1979), Bertin (1980), Lin and 
Bertin (1984), Morozov (1980), Polyachenko (1989), Griv and Peter 
(1996), Polyachenko and Polyachenko (1997), initially the disk is Jeans- 
stable against the small-scale axisymmetric (radial) perturbations but 
unstable against the relatively large-scale nonaxisymmetric (spiral) per- 
turbations. The initial vertical velocity dispersion was chosen equal to 
c z = 0.2c r . Finally, the angular velocity VLq was replaced by (Hohl, 
1971, 1972) 

The sense of disk rotation was taken to be counterclockwise, and 
units are such that G = 1 and the mass of a star m s = 1. The cutoff 
radius was chosen r cut = 0.01R. Within the framework of our model, 
dynamically young stars form the very thin disk h <C R. Therefore, at 
the start of simulations the disk thickness was chosen h = 0.05-R. A 
time t = 1 is taken to correspond to a single revolution of the initial 
disk. In all experiments the simulation was performed up to a time 
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t between 8 and 10. It should be noted here that after about three 
rotations the picture is always stabilized and practically no significant 
changes in gross properties of the model over this time are observed. 
Any difference between the results of simulations with and without 
applying the so-called quiet starts procedure to select the very regular 
initial coordinates of particles was not found. (Basically, by applying 
the method of quiet starts, one uses no random numbers in the initial 
conditions to suppress the noise level in a system. Such techniques 
have proven useful in obtaining realistic noise levels without the use 
of a large number of particles; Byers and Grewal, 1970.) Moreover, 
tests indicated that the results were insensitive to changes in other 
parameters, e.g. the number of stars in the range N = 10, 000 — 40, 000, 
the cutoff parameter in the range r cut = (0.005 — 0.05)R and the initial 
disk thickness in the range h = (0 — 0.2)R. Spiral amplitudes detected 
in our simulations are greatly in excess of an expectation from the 
level of particle noise and do not scale with N. This behavior is clearly 
inconsistent with the hypothesis that the nonaxisymmetric structure in 
the models can be attributed to collisional effects, inappropriate values 
of r cut or swing-amplified particle noise as advocated by Toomre (1990). 

Finally, slight corrections have been applied to the resultant veloc- 
ities and coordinates of the model stars so as to ensure the equilibrium 
between the centrifugal and gravitational forces, to preserve the posi- 
tion of the disk center of gravity at the origin and to include the weak 
effect of the finite thickness of the disk to the gravitational potential. 
Thus, the initial model is very near the dynamical equilibrium for all 
radii. 

The system of equations of motion (1) for N particles was in- 
tegrated by the standard Runge-Kutta method of the fourth order. 
All the particles move with the same constant Runge-Kutta time step 
At = 0.001t or b, where t or b is the orbital period. 



4. Results 

4.1. Collisions 

In order to be compared with an actual galaxy of stars, the numerical 
model of the disk should properly simulate a collisionless system. Oth- 
erwise, as it has been shown, e.g. by Morozov et al. (1985) and Griv 
et al. (2000a), a collisional disk is unstable with respect to a secular 
dissipative-type instability. The secular instability of spontaneous grav- 
ity perturbations might arise merely from the dissipation of the energy 
of regular rotation into ever larger amounts of heat when the system has 
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moderate dissipation; it results from perturbations which have negative 
energy in the dissipative medium (Morozov et ai, 1985; Fridman and 
Polyachenko, 1984, Vol. 2, p. 243). The secular instability produces 
structures completely unrelated to the effects we would like to model. 
In the spirit of the simple "molecular-kinetic" theory by Chandrasekhar 
(1960), the collisional relaxation time for a three-dimensional system 
from the point of view of deflections suffered by a test star in encounters 
with a given distribution of field objects is 



c 3 



2itG 2 mfn{ In 



(5) 



where c is the averaged velocity dispersion, mf is the "field" particle 
mass and rif its three-dimensional number density. Also lnA^ is the 
so-called Newton (or Coulomb in plasmas) logarithm, by means of 
which the long-range nature of the gravitational force is taken into 
account; to the order of magnitude In Ajy ~ lnTVf, where Nf is the total 
number of field particles. In actual galaxies lnAjv ~ 20. According 
to Eq. (5), the chosen quantities c, N, etc. guarantee that the model 
may be considered as a collisionless one to a good approximation at 
least during ~ 50 rotations. So we expect, in our numerical calculations 
because N is large, relaxation from two-body gravitational interactions 
is unimportant. We suggest that the structures observed in the A-body 
simulations originate from the collective modes of practically collision- 
less galactic models — gravitational Jeans-type modes and bending 
firehose-type modes (see below Section 4.2). 

Consider a system of mutually gravitating particles. The local 
distribution function /(r, v, t) must satisfy the Boltzmann kinetic equa- 
tion 

df df df fdf\ 



dt dr <9v \dt y coll 

where <fr(r, t) is the total gravitational potential determined self-consistently 
from the Poisson equation, 



V 2 $ = 4vrG 



J fdv, (7) 



{df/dt) co \\ oc Vcifo — f) is the so-called collisional integral which defines 
the change of / arising from ordinary interparticle collisions, v c is the 
collision frequency and fo is the quasi-steady state distribution function 
(Binney and Tremaine, 1987, p. 506; Griv et al, 2000a). 

In plasma physics, Lifshitz and Pitaevskii (1981, p. 115) have dis- 
cussed phenomena in which interparticle collisions are unimportant, 
and such a plasma is said to be collisionless (and in the lowest-order 
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approximation of the theory one can neglect the collision integral in 
the kinetic equation). It was shown that a necessary condition is that 
v c <C where oj is the typical frequency of collective oscillations: then 
the collision operator in the kinetic equation is small in comparison with 
df/dt. In the Appendix, we have shown that the frequency of Jeans 
oscillations in a stellar disk to ~ fi. Therefore, in the gravitation case 
in the lowest-order approximation of the theory we can neglect the 
effects of collisions between particles on a time scale of many rotations 
if v c <C £1. Lifshitz and Pitaevskii have pointed out that collisions may 
be neglected also if the particle mean free path is large compared with 
the wavelength of collective oscillations. Then the collision integral in 
Eq. (6) is small in comparison with the term v • {df/dr). 

In this subsection, we test numerically if the models used in our 
iV-body simulations are being correctly modeled as collisionless Boltz- 
mann (Vlasov) systems. The direct method of checking if the system 
is being modeled as a collisionless system is to repeat a calculation 
using a mass spectrum (Rybicki, 1971). It is obvious that as a result of 
gravitational collisions there is a tendency towards energy equipartition 
between the various masses. Hohl (1973) has determined the experi- 
mental relaxation time and compared it with a theoretical prediction 
for the collisional relaxation time of a two-dimensional disk. Here we 
do such a comparison by using the three-dimensional system. 

Let us consider the three-dimensional computer model consisting 
of 2% stars of mass 7713 = 10m s , 18% stars of mass 777-2 = 2?77 s and 
80% stars of mass m\ = 0.55m s . The total number of stars, which are 
distributed in the system in accord with the law (2) is N = 20, 000. 
Initially, the different mass groups of stars are distributed with the 
same velocity dispersion (with different temperatures). 

In the Appendix of the present paper we demonstrate that the 
disk becomes almost stable gravitationally for c r > 2ct- In such a 
Jeans-stable system collective effects associated with the classical grav- 
itational instability will not affect the random velocity dispersion of 
particles: then the change of velocity dispersion can be explained only 
by usual two-body encounters. For this reason the initial condition was 
chosen to be a quasi-stable uniformly rotating disk with c r = 2ct- 
Following Chandrasekhar (1960) and Hohl (1973), let us define the 
relaxation time te as the time required for the mean change of the 
kinetic energy per unit mass of the test star to equal the initial kinetic 
energy. 

In Figure 1 we show the change of the ratio of the mean particle 
(kinetic) energy, < m\Vl >, < m^V^ 2 > and < TO3V3 2 > (in units of 
the total kinetic energy of the system), for the different mass groups, 
where Vi is the total velocity of a given mass group. As is expected, the 
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Figure 1. Rate of change of the mean kinetic energy for stars of the three mass 
groups of the Jeans-stable model with N = 20, 000 and c r = 2ct; 1 — kinetic energy 
of stars with the mass of a star 10m s , 2 — with the mass of a star 2m s and 3 — with 
the mass of a star 0.55m s . The two groups of heavy stars lose kinetic energy while 
the group of lightest stars gains an approximately corresponding amount of kinetic 
energy. In general agreement with the analytical estimation, the mean slope of the 
curves will result in energy equipartition after about 100 — 200 rotation periods. 
This result shows that the model used in our iV-body simulations is "collisionless" 
to a good approximation and suggests that interparticle collisions do not play a 
significant role for instabilities studied in the paper. 



two groups of heavy stars lose energy while the group of lightest stars 
gains an approximately corresponding amount of kinetic energy Also 
as is expected, one can see the decrease in the change of the kinetic 
energy with time. This is because the collisional frequency v c « 1/te is 
inversely proportional to the velocity dispersion (Eq. [5]), and thus the 
encounters only weakly affect the stars with high random velocities. 

As one can see, the mean slope of the curves shown in Figure 
1 will result in energy equipartition after about 100 — 200 rotation 
periods. It is crucial to realize that these relaxation times even for 
this relatively small number of model stars are much longer that the 
time of a single disk revolution. We conclude that the computer models 
used in the present study may be considered as collisionless ones to a 
good approximation at least during the first 10 rotations which are of 
especial interest in spiral-galaxy simulation. Therefore, we argue that 
the collective effects studied in the paper were apparent before the 
collisional time scale associated with the binary interactions of stars 
was reached. 



4.2. Warm Toomre-stable disk 

To study the development of spiral structure and to confirm the theo- 
retical predictions made in the Appendix, warm disks with initial radial 
velocity dispersion of identical stars equal to Toomre's (1964) stabiliz- 
ing velocity dispersion ct were modeled. Figure 2 shows an example 
of the development of spiral structures. The structure grows rapidly 
on a dynamical time scale ~ f2 _1 . At a time t ~ 0.6 a multi-armed 
filamentary spiral structure forms with several short, tightly wound 
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Figure 2a. The time evolution of a warm Toomre-stable three-dimensional disk of 
stars. The system is violently unstable to low-m spiral Jeans-type modes. After ~ 3 
rotations the spiral structure almost dissapears, and a prominent bar forms. 



Figure 2b. Continued 



trailing arms.[] Soon after, at t ~ 2.8 a bar structure appears with 
two almost symmetrical open spiral arms. However, sometimes, e.g. at 
t = 1.6 or at t = 2.6, one can see a "one-armed" galaxy consisting of 
a prominent arm and a short, diffuse second one. After ~ 3 rotations, 
the spiral arms almost disappear, and a bar forms in the central parts 
of the system. 

The m = 1 instability shifts the point with highest density from 
the center of mass. It is interesting to note that in many disk-shaped 
galaxies, e.g. in the spiral galaxies M 101 and NGC 1300, there ap- 
pears to be a deviation from rotational symmetry. In principle, such a 
deviation may be due to this one arm instability. Note that the latter 
has already been suggested by B. Lindblad (see Toomre, 1996). In fact, 
asymmetries in the distribution of light in spiral galaxies have been 
known for a long time. Baldwin et al. (1980) presented a number of 
spiral galaxies in which the distribution of neutral gas is lopsided, i.e. 
the gas extends further out on one side of the galaxy than on the other. 
In many systems asymmetries are also found to be well developed in the 
stellar disks as indeed confirmed recently by near-infrared observations 
(Zaritsky and Rix, 1997). Rudnick and Rix (1998) quantified the mean 
asymmetry of 54 early-type disk galaxies using the amplitude of the 
m = 1 azimuthal Fourier component of the i?-band surface brightness. 
It was found that 20% of all disk galaxies have azimuthal asymmetries 
in their stellar disk disribution, confirming lopsidedness as a dynamical 
phenomenon. Based on the iV-body simulations, we can offer a tentative 
hypothesis which accounts lopsided asymmetries in spiral galaxies by 
the m = 1 Jeans-unstable spiral mode. 

Following Athanassoula and Sellwood (1986), in Figure 3 the time 
evolution of the Fourier spectrum of the spiral pattern shown in Figure 

1 Interestingly, according to fairly recent observations, such ragged galaxies with 
almost chaotic-looking, "flocculent" spiral structures are more common than the 
grand design ones (Binney and Tremaine, 1987, p. 391). 
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Figure 3. The time evolution of the Fourier spectrum of the particle distribution 
shown in Figure 2 for different azimuthal mode numbes m. In this figure, the ordinate 
is the Fourier coefficient A(p, m) and the abscissa is the pitch angle p. Spiral ampli- 
tudes are greatly in excess of an expectation from the level of particle noise 
true instabilities develop in the model. The result suggests that very flat stellar disks 
of galaxies can be unstable to the most unstable one-armed (m = 1), two-armed 
(m = 2) and three-armed (m = 3) spiral Jeans-type gravity perturbations. The 
one-armed mode may explain the asymmetric patterns observed in many isolated 
galaxies. 

2 is plotted. Only the azimuthal m = 1 — 6 components are shown. 
The complex Fourier coefficients, A, are determined from the following 
summation 

1 N 

A(p,m) = — ^exp{i[m<^- +pln(r,-)]}, 

where N is the number of particles and (ipj,rj) are the polar coordi- 
nates of the j-th particle. The pitch angle of an m-armed logarithmic 
spiral ip is given by ip = arctan(p/m), and positive p corresponds to 
trailing spirals and negative p to leading ones. 

This logarithmic spiral Fourier analysis shows that the peaks of 
the signals move to positive p with increasing time, reach maxima at 
t = 1.0 — 1.4, and eventually decay. It follows, then, that after a time 
t = 2.0 — 3.0 the model becomes practically stable with respect to 
gravitational Jeans-type modes. We can observe that the m = 1 — 3 
modes are the most unstable ones. An interesting feature is that all 
these modes have roughly the same pitch angle. The amplitudes of the 
signals grow aperiodically, supporting the assumption that we have deal 
with the aperiodic Jeans- type instability of gravity disturbances. 

In actual galaxies, discrete spiral modes were already found by Rix 
and Zaritsky (1995), Block and Puerari (1999) and Block et al. (1999). 
It was stated that in the near-infrared, the morphology of older star- 
dominated disk indicates a simple classification scheme: the dominant 
Fourier m-mode in the dust penetrated regime, and the associated pitch 
angle. A ubiquity of low m = 1 and m = 2 modes was confirmed. 
Fridman et al. (1998) detected m = 1 — 8 spiral modes in relatively 
young stellar population of the nearly face-on galaxy NGC 3631 from 
observations in the H Q line. 

From the edge-on view pictured in Figure 4, one can see that 
from an initial very thin model, a fully three-dimensional disk develops 
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Figure 4a. Higher-resolution plots of the central parts (edge-on view) for the simu- 
lation run seen in Figure 2. At a time t « 1.4, one can see the rapid development of 
the snake-shaped, bending firehose-type instability as a precursor of galactic bulge 
formation in the central, rigidly rotating portions of the system. This instability is a 
transient one, and develops on a time scale of single rotation of the system only. The 
linear dimension of the bending along the z-axis is comparable to the disk thickness 
and much less than the radius of the outermost parts of the system. The usual name 
for this instability of collective systems — firehose instability — recalls the fact that 
the unstable motion of a hose when the flow velocity of water inside becomes very 
high has essentially the same nature. 



Figure 4b. Continued 



immediately at t pa 0.4 with a mean height < z > above the plane, 
corresponding to the force balance between the gravitational attraction 
in the plane and the "pressure" due to the velocity dispersion (i.e. 
"temperature" ) in the vertical, z-direction. A straightforward estimate 
shows that 

where (i is the three-dimensional mass density in the mid-plane. After 
a time t pa 0.4 there is no change in the edge-on structure until at 
t pa 1.4 (when the Jeans instability is almost switched off by increasing 
the planar velocity dispersion). At somewhat later times, t pa 3.2, 
a "box-shaped" or "peanut-shaped" bar structure is developed. The 
development of the box/peanut bar structure is explained in the terms 
of the bending firehose-type instability (Raha et al., 1991; Merritt and 
Sellwood, 1994; Griv and Chiueh, 1998). 

One sees that at times t = 1.4 — 2.4, the bending instability fiercely 
develops in the central, almost rigidly rotating parts of the system and 
is switched off at t pa 2.8. At later times, no dramatic evolution is 
observed in our simulations. To emphasize, this type of bending in- 
stabilities develops in the central, rigidly rotating parts of the system. 
The three-dimensional simulations show that the bending instabilities 
mentioned above destroy the strong bar formation (Griv and Chiueh, 
1998). Note that the bending modes have been found experimentally 
by Merritt and Hernquist (1991) and Merritt and Sellwood (1994), 
but in the nonrotating models. Combes et al. (1990) and Raha et al. 
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(1991) found the bending instability in the central regions of a fast- 
rotating iV-body disk. It was stated that this instability may play a 
part in the formation of triaxial bulges in flat galaxies and explanation 
was given why many bulge stars are less than 10 Gyr old (Raha et 
ai, 1991). Bureau and Freeman (1997) have argued that new observa- 
tions constitute a strong case in favor of the bending (or bar-buckling) 
mechanism for the formation of boxy/peanut-shaped bulges in spiral 
galaxies. Observations suggest very strongly that many of ordinary 
spiral galaxies, including the Milky Way, have similar small triaxial 
bars in their central regions (Dwek et al, 1995). 

Apparently, this instability of a sufficiently thin stellar disk has 
been predicted by Toomre (1966) by using the simplified theory based 
on moment equations. Toomre considered the collisionless analog of the 
Kelvin-Helmholtz instability in an infinite, two-dimensional, nonrotat- 
ing sheet of stars. This is the usual way to discuss the conditions of 
the firehose instability in plasma physics (Krall and Trivelpiece, 1986). 
It was demonstrated by Toomre that the bending instability is driven 
by the stellar "pressure" anisotropy: the source of free energy in the 
instability is the intrinsic anisotropy of a velocity dispersion ("tem- 
perature"). The bending firehose- type instability was also discovered 
independently by Kulsrud et al. (1971) and Mark (1971) with a more 
accurate kinetic theory. Kulsrud et al. have described the mechanism 
of the instability in terms of the balance of centrifugal force acting on 
the stars of a bending layer and the gravitational attraction toward the 
plane of the practically nonrotating, pressure-supported disk. Fridman 
and Polyachenko (1984) have discussed the role of the instability in 
explaining the existence of maximum oblateness in almost nonrotat- 
ing elliptical galaxies and the formation of the bulges of disk-shaped 
galaxies of stars. 

As is seen in Figures 2 and 4, the basic form of warm disk evolution 
is the expansion of the outermost parts together with the contraction 
of a large part of the mass towards the center; in the final quasi-steady 
state the surface density decreases exponentially outwards (Fig. 5). In a 
future publication we intend to explain the exponential surface density 
distribution of the final quasi-steady state model by applying the quasi- 
linear theory of a stellar disk's stability. Note only that the study of 
surface photometry have been shown that most spiral and SO galaxies 
have an exponential disk with radial surface-brightness distribution ~ 
exp(— const • r) (Freeman, 1970). 

At first, the surface density and rotation velocity change rapidly 
(Figs. 5 and 6), but after a time t ~ 1.0 there is little further change 
during the following two rotations, except in the central portions of the 
disk. The model generates a nonflat rotation curve (Fig. 6); in addition 
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Figure 5. Variation of the surface density for a disk shown in Figure 2. Note the 
rapid change with time of the surface density at t < 1. The basic form of the disk 
evolution is the expansion of the outermost portions together with the contraction 
of a large part of the mass towards the center. In a quasi-steady state after about 
2 — 3 rotations, the surface density decreases exponentially outwards. According to 
observations, stellar disks of most spiral and SO galaxies are similar approximated 
well by a function that is exponential in r. 



Figure 6. The azimuthal average rotation curve for a stellar disk shown in Figure 
2. Note the non-flat rotation curve at the periphery of the system; growing density 
waves (spiral arms) transfer the angular momentum of the circular motion outwards. 



growing density waves (spiral arms) transfer the angular momentum of 
the circular motion outwards (see the difference in the rotation curves 
at t = 1.0 and t = 2.0 in Figure 6). As mentioned above, the model 
practically shows no the time evolution of the surface density and the 
velocity curve after t ~ 2. The latter seem to indicate that at that time 
the system has nearly reached a quasi-steady state (cf. Hohl, 1971). 

The time evolution of the velocity dispersions c r , c v and c z is shown 
in Figure 7. 

Initially, during the first rotation the velocity dispersions along 
each axis are changed rapidly. During the following two rotations the 
velocity dispersions are changed with time slightly, confirming that 
between t = 2 and t = 3 the system has nearly reached a quasi-steady 
state. However, in the central portions of the disk the velocity dispersion 



Figure 7. Variation of velocity dispersions for a disk shown in Figure 2; 1 — radial 
velocity dispersion, 2 — azimuthal velocity dispersion and 3 — vertical velocity 
dispersion. Initially, during the first rotation the velocity dispersions increase with 
time rapidly; later the dispersions grow only slightly. However, in the central portions 
of the disk the vertical velocity dispersion c z grows rapidly even after the time 
t ~ 1.5, when the bending firehose-type instability begins to to operate (see Figure 
4). In this central region the bending instability strongly changes the ratio of vertical 
to radial (and tangential) dispersion by transfering kinetic energy of plane random 
motions of particles to vertical motions. 
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Figure 8. Comparison of Toomre's radial-velocity dispersion and the modified crit- 
ical dispersion with the experimental radial velocity dispersion of a disk shown in 
Figure 2 at t = 2.5; f — the experimental radial velocity dispersion, 2 — the 
modified velocity dispersion, 3 — Toomre's velocity dispersion. 



Figure 9. The time development of the mean-square spread in planar random 
velocity v\ of model particles for different radii r. 



c 2 increases rapidly even after the time t « 1.5, when the bending 
instability begins to operate. 

As was stated, one of the goals of the present research is to com- 
pare Toomre's (1964) critical dispersion (4) and the modified dispersion 
given by Eq. (25) with the experimental radial velocity dispersion of 
the quasi-steady state model of a stellar disk. In what follows, we are 
doing such a comparison. 

As is seen in Figure 8, the modified critical dispersion, in contrast 
to Toomre's (1964) critical dispersion, is in fair agreement with the 
experimental data, especially for the outer, differentially rotating parts 
of the disk. The agreement is quite reasonable considering the rough 
nature of the theory. In the central parts of the model the agreement 
with the theory is not good. This is to be expected since the condition 
of nearly circular orbits adopted in the Appendix does not fulfill in the 
central almost nonrotating regions of model disks (and in spiral galax- 
ies, including both barred and normal galaxies, and in disks of elliptical 
galaxies). Therefore in the central parts of the disk (r < 0.4 — 0.5) 
the condition for the local criteria (both the Toomre and the modified 
criterion) to be applicable is violated. 

One of the important problems of stellar disks is the determination 
of the random velocity diffusion. Such velocity diffusion can be caused 
by stellar disk turbulence. Velocity diffusion produced by turbulence 
tells us something about the turbulence process. 

To compute velocity diffusion we calculate the mean-square spread 
in the planar random (residual) velocity v \ = + (? as a function of 
time of model particles for different radii r. In general v\ has a time 
dependence like that shown in Figure 9. 

As is seen, in the differentially rotating (r > 0.4 — 0.5), Jeans- 
unstable (at times t < 1.8 — 2.0) parts of the disk approximately v\ oc 
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Figure 10. The time evolution of a hot three-dimensional disk of stars. The system 
is stable to all Jeans modes. The result agrees with the theoretical explanation 
described in the paper. 

t. This behavior has been predicted analytically (see the Appendix). 
Interestingly, observations indicate a similar increase of v\ with time 
in the solar neighborhood (Binney and Tremaine, 1987, p. 470; Gilmore 
et al, 1990). 

4.3. Hot Jeans-stable disk 

The time evolution of the hot disk was investigated, in which the initial 
dispersion of random radial velocities of model stars was chosen to be 
greater than Toomre's critical dispersion, namely, c r = 2ct (see Eq. 
[24]). For such calculations we obtain plots like those shown in Figure 
10. In agreement with the theory, the system becomes gravitationally 
(Jeans-) stable. 



The linear Lin-Shu density wave theory is considered now as well es- 
tablished (Lin and Shu, 1966; Lin et al, 1969; Shu, 1970; Bertin, 1980; 
Binney and Tremaine, 1987, p. 347; Griv and Peter, 1996; Griv et al, 
2000b). However, as soon as the amplitude of the gravity perturbations 
becomes appreciable a new generation of problems, namely, nonlinear 
problems, become important. Below we present a quantitative weakly 
nonlinear theory of the classical Jeans instability of small-amplitude 
gravity disturbances (e.g. those produced by a barlike structure, a spon- 
taneous perturbation and/or a companion galaxy) in self-gravitating, 
rapidly rotating stellar disks of flat galaxies. 

In the rotating frame of a disk galaxy, the collisionless motion 
of an ensemble of identical stars in the plane of the system can be 
described by the Boltzmann equation for the distribution function 
f(r,v,t) without the integral of collisions (Lin et al, 1969): 



5. Appendices 




0, (9) 



grivetalpaper.tex; 1/02/2008; 16:07; p. 17 



18 



E. GRIV ET AL. 



where r, (p, z are the galactocentric cylindrical coordinates, the total 
azimuthal velocity of the stars was represented as a sum of the random 
and the basic rotation velocity V ro t = rf2, and v r is the random 
velocity in the radial direction. The quantity £l(r) denotes the angular 
velocity of galactic rotation at the distance r from the center, and 
the epicyclic frequency k varies from 2f2 for a rigid rotation to fi 
for a Keplerian one. Random velocities are small compared with r£l. 
Collisions are neglected here because the collision frequency is much 
smaller than the cyclic frequency. In the kinetic equation (9), $(r,i) is 
the total gravitational potential determined self-consistently from the 
Poisson equation (7). 

The equilibrium state is assumed, for simplicity, to be an ax- 
isymmetric and spatially homogeneous stellar disk. Secondly, in our 
simplified model, the perturbation is propagating in the plane of the 
disk. This approximation of an infinitesimally thin disk is a valid ap- 
proximation if one considers perturbations with a radial wavelength 
that is greater than the typical disk thickness. We assume here that 
the stars move in the disk plane so that v z = 0. This allows us to use 
the two-dimensional distribution function / = f(v r ,v v , t)S(z) such that 
J fdvrdv^dz = a, where c is the surface density. We expect that the 
waves propagating in the disk plane have the greatest influence on the 
development of structures in the disk. The latter suggestion is strongly 
supported by numerical simulations (Hohl, 1978). 

The disk in the equlibrium is described by the following equation: 



where df e /dt = and the angular velocity of rotation £l(r) is such that 
the necessary centrifugal acceleration is exactly provided by the central 
gravitational force r$7 2 = d& c /dr, where $ e is the axisymmetric equi- 
librium potential. Equation (10) does not determine the equilibrium 
distribution f c uniquely. For the present analysis we choose / e in the 
form of the anisotropic Maxwellian (Schwartzschild) distribution 



The Schwarzschild distribution function is a function of the two epicycli 
constants of motion £ = v\/2 and rgfi(ro), where r$ = r + (2Q/n 2 )v< fi 
These constants of motion are related to the unperturbed star orbits: 

r = [sin(0o — nt) — sin </>q] ; v r = v± cos (<po — nt) ; 




(10) 




(11) 



K 
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2Q v± r . , 

if = COS((po — Kt) — COS 00 ; 

k r K 

dip v_\_d£l k . . , . , , 



where v±, 4>o are constants of integration, v±/ nr$ ~ p/fQ <C 1, p is the 
mean epicycle radius, and we follow Lin et al. (1969), Shu (1970) and 
Lynden-Bell and Kalnajs (1972), making use of expressions for the un- 
perturbed epicyclic trajectories of stars (12) in the equilibrium central 
field < & e (r). In Eqs. (12), tq is the radius of the circular orbit, which is 
chosen so that the constant of areas for this circular orbit r^dipo/dt) 
is equal to the angular momentum imtegral M z = r 2 (dip/dt), and 
v 2 ^ = v 2 + (2^1/k) 2 v' 2 . Also, ipo is the position angle on the circular orbit, 
(dtpo/dt) 2 = (l/ro)(d$ e /dr)o = £l 2 . The quantities k and c r are eval- 
uated at ro- In Eq. (11) the fact is used that as follows from Eqs. (12) in 
a rotating frame the radial velocity dispersion c r and the azimuthal ve- 
locity dispersion c<o are connected through c r ~ (2CI/k)cu,. In the Solar 
vicinity, 20, /k ~ 1.7. The distribution function f e has been normalized 
according to f^f^fedvrdvy, = 2-k(k/2Q) / °° v ± dv±f c = a e , where 
a c is the equilibrium surface density. Such a distribution function for 
the unperturbed system is particularly important because it provides 
a fit to observations (Lin et al, 1969; Shu, 1970). 

We proceed by applying the standard procedure of the quasi-linear 
(weakly nonlinear) approach (Lifshitz and Pitaevskii, 1981; Alexandrov 
et al., 1984; Krall and Trivelpiece, 1986) and decompose the time depen- 
dent distribution function / = /o(v, t) + /i(v, t) and the gravitational 
potential $ = $o(r,t) + $i(r,i) with |/i// | < 1 and |$i/$ | < 1 
for all r and t. The functions f\ and $i are oscillating rapidly in 
space and time, while the functions fo and <E>o describe the slowly 
developing "background" against which small perturbations develop; 
fo(t = 0) = f e and ^>o(^ = 0) = <3? e - The distribution /o continues to 
distort as long as the distribution is unstable. Linearizing Eq. (9) and 
separating fast and slow varying variables one obtains: 

dh = &hdfo + ld$idfo (13) 

dt dr dv r r dip dv^ 

df _ d^dh | l^ifl/i, (14) 
dt dr dv r r dip dv v ' 

where d/dt means the total derivative along the star orbit (12) and 
(• • •) denotes the time average over the fast oscillations. To emphasize it 
again, we are concerned with the growth or decay of small perturbations 
from an equilibrium state. 
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In the epicyclic approximation, the partial derivatives in Eqs. (13) 
and (14) transform as follows (Lin et al, 1969; Shu, 1970; Morozov, 
1980; Griv and Peter, 1996): 

d d 20. v v d d __ (20\ 2 d 20 v r d 

dv r r d£ k v\ d(f>o ' dv v \ k J ^88 k v \ dfo ' 

To determine oscillation spectra, let us consider the stability prob- 
lem in the lowest WKB approximation; this is accurate for short wave 
perturbations only, but qualitatively correctly even for perturbations 
with a longer wavelength, of the order of the disk radius R. In this 
local WKB approximation in Eqs. (13) and (14), assuming the weakly 
inhomogeneous disk, the perturbation is selected in the form of a plane 
wave (in the rotating frame): 

where 8f, 5& are amplitudes that are constant in space and time, m 
is the nonnegative azimuthal mode number (= number of spiral arms), 
lu* = oj — mfl is the Doppler-shifted wavefrequency and \k r \R S> 1 (Lin 
et al, 1969; Shu, 1970; Griv and Peter, 1996). The solution in such a 
form represents a spiral wave with m arms whose shape in the plane is 
determined by the relation 

K{r ~ r ) = -m((p - tp ). (17) 

With ip increasing in the rotation direction, we have k r > for trailing 
spiral patters, which are the most frequently observed among spiral 
galaxies. A change of the sign of k r corresponds to changing the sense 
of winding of the spirals, i.e. leading ones. With m = 0, we have the 
density waves in the form of concentric rings that propagate away from 
the center when k r > or toward the center when k r < 0. 

In Eq. (13) using the transformation of the derivatives d/dv r and 
d/dvy given by Eqs. (15), one obtains the solution 

/l = .L^aFaT' (18) 

where f\{t' = — oo) — > 0. In this equation making use of the time 
dependence of perturbations in the form of Eq. (16) and the unper- 
turbed trajectories of stars (12) in the exponential factor occurring in 
the formula (18), it is straightforward to show that 

flf n ~ p j(™-0W>o-C) j,(y)J (y) 

h = -*iwf E E <* £ „_ J ;1 X)J ° (X) . (w) 

i=-oon=-oo 
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where J/(x) is the Bessel function of the first kind of order I, x = 
k*v±/K ~ k*p, Zc* = k{l + [(20/k) 2 - 1] sin 2 tp} 1 / 2 is the effective 
wavenumber, ip is the pitch angle of perturbations, tan^ = k v /k r = 
m/rk r , and we used the expansion 



and the usual Bessel function recursion relation J/ + i(%) + Ji-\{x) = 
(2l/x)Ji(x)- in Eq. (19), the denominator vanishes when — Ik — > 
0. This occurs near corotation (Z = 0) and other resonances (Z = 
±1,±2, •■■). The Lindblad resonances occur at radii where the field 
(<9/(?r)<3?i resonates with the harmonics I = — 1 (inner resonance) and 
Z = 1 (outer resonance) of the epicyclic (radial) frequency of equilibrium 
oscillations of stars k. Clearly, the location of these resonances depends 
on the rotation curve and the spiral pattern speed; the higher the m 
value, the closer in radius the resonances are located (Lin et al, 1969; 
Shu, 1970). In this paper, only the main part of the galactic disk is 
studied which lies sufficiently far from the resonances: below in all 
equations lu* — Ik / 0. 

In the WKB approximation, the linearized Poisson equation (7) 
is readily solved to give the improved Lin-Shu expression for the per- 
turbed surface density of the two-dimensional disk (Lau and Bertin, 
1978; Lin and Lau, 1979): 



On the other hand, integrating Eq. (19) over velocity space f dv r f dv^fi = 
7t(k/2Q) J °° fidv 2 ^ = cti, and equating the result to the perturbed den- 
sity given by the improved Lin-Shu solution (Eq. [20]), the dispersion 
relation oj* = ^(fc*) is obtained: 



where |cj*| / Ik, h{x) is a Bessel function of imaginary argument 
with its argument x = k^c^/K 2 ~ k 2 p 2 , This dispersion relation gen- 
eralizes the Lin-Shu-Kalnajs one for nonaxisymmetric perturbations 
(tp / 0) propagating in a homogeneous disk excluding the resonance 
zones (a>* — Ik ^ 0). Analyzing Eq. (21), it is useful to distinguish 
between two limiting cases: the cases of of epicyclic radius that is small 
compared with wavelength, x < 1, and of epicyclic radius that is large 
compared with wavelength, x S> 1. We study the problem when the 
frequency of disk oscillations is smaller than the epicyclic frequency. 



(X) 




n=— oo 



ai = -|A;|$i/27rG. 



(20) 




(21) 
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In the opposite case of the high perturbation frequencies, > k, 
the effect of the disk rotation is negligible and therefore not relevant 
to us. This is because in this "rotationless" case the star motion is 
approximately rectilinear on the time and length scales of interest which 
are the wave growth/damping periods and wavelengths, respectively 
(cf. Alexandrov et al, 1984, p. 113). Thus, the terms in series (19) 
and (21) for which |/| > 2 may be neglected, and consideration will 
be limited to the transparency region between the turning points in a 
disk, i.e. between the inner and outer Lindblad resonances. The I = 
harmonic that defines the corotation resonance dominates in series. 

Thus, |a;*| less than the epicyclic frequency of any disk star. In this 
case, in Eq. (21) function A(ar) = exp(— x)Ii(x) starts from A(0) = 0, 
reaches a maximum A max < 1 at x ~ 0.5, and then decreases. Hence, 
the growth rate of the instability has a maximum at x < 1 (Griv et 
al., 1999b, Fig. 2 in their paper). In addition, we consider only low- 
m perturbations (which are important only for the problem of spiral 
structure): from now on in all equations m ~ 1. 

Using the well-known asymptotic expansions of Bessel functions, 
the simplified (|/| < 1) dispersion relation reads 



where Wj « k 2 - 2TrGao\k\F(x) is the squared Jeans frequency and 
F « Ik? e~ x I\{x) I k 2 <? r is the so-called reduction factor, which takes 
into account the fact that the wave field only weakly affects the stars 
with high random velocities. The latter means that the spiral structure 
is seen more distinctly in objects with the lowest velocity dispersion. In 
the limit x < 1, F(x) « (A;*/A;) 2 [1 — x + (3/4)x 2 ] (but, of course, in order 
to be appropriate for a WKB wave we consider the perturbations with 
\k r \r > 1). In the opposite limit x > 1, F(x) ps (l/kp) 2 [l-(l/2irx) 1 / 2 ]. 
Different forms of F(x) are given by Athanassoula (1984). The aperi- 
odic Jeans instability (gravitational collapse) is possible when u) 2 = 
ijj 2 < 0. That is, the real part of the frequency of this unstable motion 
vanishes in a rotating frame. The Jeans instability is driven by a strong 
nonresonant interaction of the gravity fluctuations with the bulk of the 
particle population, and the dynamics of Jeans perturbations can be 
characterized as a fluid-like interaction. Therefore, such an instability, 
which implies the displacement of macroscopic portions of a gravitating 
system, may be also analyzed through the use of relatively simple hy- 
drodynamic equations (Lau and Bertin, 1978; Lin and Lau, 1979; Lin 
and Bertin, 1984). In general, the growth rate of the Jeans instability is 
relatively high |3u>j| ~ Q; perturbations with wavelength Aj « 2ttc t /k 
have the fastest growth rate (Griv et al, 1999a, b; Griv et al, 2000a, b). 
In the Solar vicinity of the Galaxy, Aj = 2 — 3 kpc. 




(22) 
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Equation (22) has two roots which describe two branches of oscil- 
lations: long-wavelength, x < 1, and short- wavelength, x S> 1, Jeans 
branches. The dispersion law of low-m Jeans oscillations is 

<^*i,2 = ±p|a;j|. (23) 

Here p = 1 for Jeans-stable (u; 2 > 0) perturbations and p = i for 
Jeans-unstable (w 2 < 0) perturbations. 

The Jeans perturbations can be stabilized by the random velocity- 
spread. Indeed, if one recalls that such unstable perturbations are pos- 
sible only when oj^ 2 = w 2 < 0, then by using the condition ujj > for 
all possible k, from Eq. (22) a local stability criterion to suppress the 
instability of arbitrary Jeans-type gravity perturbations can be easily 
obtained (Griv et al, 1999a,b; Griv et al, 2000b): 

c r > c r ,crit = c T {l + (2^/k) 2 - 1] sin 2 V} 1/2 , (24) 

where c r is the radial- velocity dispersion, c r)Cr it is the critical radial 
velocity dispersion and usually in the main part of a disk-shaped galaxy 
the parameter 2Q/k = 1.5— 1.7. It is clear from this local criterion that 
stability of nonaxisymmetric Jeans perturbations (sin tp ^ 0) in a differ- 
entially rotating disk (2Q/k > 1) requires a larger velocity dispersion 
than Toomre's critical value op. It is only for rigidly rotating disks 
(20,/ k = 1) and/or for axisymmetric perturbations (sin?/> = 0) that 
the critical value of c r is ct- Hence, the ordinary Toomre's velocity 
dispersion ct should stabilize only small-scale axisymmetric (radial) 
perturbations of the Jeans type. The differentially rotating disk however 
is still unstable against relatively large-scale nonaxisymmetric (spiral) 
modes with -0 7^ 0, in particular very open modes with ip > 45°. Accord- 
ing to Polyachenko (1989) and Polyachenko and Polyachenko (1997) 
the marginal stability condition for Jeans perturbations of an arbitrary 
degree of axial asymmetry has been available since 1965 (Goldreich and 
Lynden-Bell, 1965), though in a slightly masked form. A relationship 
exists between Eq. (24) and what Toomre (Toomre, 1981; Binney and 
Tremaine, 1987, p. 375) called "swing amplification." 

Summarizing, for Jeans-stable differentially rotating stellar disks, 
the critical value of the radial velocity dispersion must be greater than 
(although of the order of) ct- In accordance with Eq. (24), 

c M ~ {2Q/k)c t ~ 2c T (25) 

might be the "new" Toomre's velocity dispersion to suppress all Jeans 
perturbations in a stellar disk, including the most dangerous open ones. 
Note that the destabilizing effect of tangential (ip 7^ 0) gravitational 
forces (pitch angle dependent effect) was first considered by Lau and 
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Figure 11. The generalized Lin-Shu-Kalnajs dispersion relation of a homogeneous 
stellar disk in the case 20, /k = \/2 and | sini/>| = 1 for the different radial velocity 
dispersion values: (a) c r = 0.5cm, (b) c r = 0.8cm, (c) c r = cm and (d) c r = 1.5cm- 
The solid curves represent the real part of the dimensionless Doppler-shifted 
wavefrequency v = uj*/k of low- frequency (|ai*| < n) long- wavelength (1) and 
short-wavelength (2) Jeans oscillations we are interested in. The dashed curves 
represent the imaginary part of the dimensionless wavefrequency of low-frequency 
vibrations. The dot-dashed curves represent the wavefrequencies of additional 
high-frequency > k) Jeans modes. 



Bertin (1978) and Lin and Lau (1979) using a much simpler gaseous 
approach, and a stability criterion that is similar to Eq. (25) was also 
derived. Clearly, the criterion (25) is only approximate one. In the 
present study, we test the validities of this modified criterion for sta- 
bility of Jeans perturbations in a self-gravitating, infinitesimally thin 
and almost collisionless disk of stars. 

Correspondingly, for Jeans-stable differentially rotating stellar disks, 
the critical value of the widely used Toomre's instability parameter 
Q = c rjCr it/cT must be greater than 20/k ~ 2. (Toomre's Q- value 
is a measure of the ratio of thermal and rotational stabilization to 
self-gravitation.) Observations already demonstated that stellar (and 
gaseous) disks of a large number of galaxies, including Milky Way's disk, 
are near the boundary of the gravitational stability with Q between 
2 and l\ over a large range of galactic disks (Fridman et ai, 1991; 
Bottema, 1993). 

A general impression of how the spectrum of nonaxisymmetric 
Jeans perturbations behaves in a homogeneous nonuniformly rotating 
disk can be gained from Figure 11, which shows the dispersion curves 
in the cases of Jeans- unstable systems ((a) and (b)), a marginally 
Jeans-stable system (c) and a Jeans-stable one (d) for values of I = 

0, ±1, ±2, ±3 (as determined on a computer from Eq. [21]). 

In this figure, the ordinate is the effective wavenumber k* measured 
in terms of the inverse epicyclic radius p and the abscissa is v = uj*/k, 

1. e. the dimensionless angular frequency at which the stars meet with 
the pattern, measured in terms of the epicyclic frequency k. In general, 
for fixed wavefrequency v there are two solutions in k*p, comprising a 
long- wavelength wave, k*p < 1, and a short- wavelength wave, fc*p > 
1. A property of the solution (21) is that in a homogeneous system 
the Jeans-stable modes those with c r > cm are separated from each 
other by frequency intervals where there is no wave propagation: gaps 
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Figure 12. The growth rate of the Jeans instability of a homogeneous stellar disk 
(arbitrary units) in the case 20./ k = \/2, ro = 8, \k r \ — 1, p = 1 and | sin t/> = 1 
for the different values of the azimuthal mode number m. The growth rates have 
maxima with respect to the critical mode number m cr it ~ 1. The low-m spiral modes 
(m < 10 and m ^ 0) are more unstable than the radial ones (m = 0) and the high-m 
ones (m > 10). 



occur between each harmonic (cf. the Bernstein modes in a magnetized 
plasma) . 

The growth rate of the Jeans instability of a homogeneous stellar 
disk as determined on a computer from the expression y/2irGao\k\F(x) 
is shown in Figure 12. 

As one can see visually in this Figure, the growth rates have max- 
ima with respect to the critical mode number m cr it ~ 1. It means 
that of all harmonics of initial perturbation, one perturbation with 
the maximum of the growth rate and with m = m cr jt will be formed 
asymptotically in time. The low-m spiral modes (m < 10 and 
are more unstable than the radial ones (m = 0) and the high-m ones 
(m > 10). These low-m spiral modes are only important in the problem 
of galactic spiral structure because in contrast to the high-m modes, 
they do extend essentially over a large range of the galactic disk (Lin 
et al, 1969; Shu, 1970; Griv and Peter, 1996). 

The shape and the number of spiral arms depend on the equilib- 
rium parameters of a galaxy. For the Galaxy the most probable pattern 
is that of 1 — 4 arms, the radial distance between the arms being 2 — 4 
kpc. Such a pattern corresponds to a spiral wave mode with maximum 
of growth rate. 

Next we substitute the solution (19) into Eq. (14). Taking into 
account that the terms I / n vanish for axially symmetric functions /o, 
after averaging over 4>q we obtain the equation for the slow part of the 
distribution function: 

§ = »EEi^i 2 7^— £2^. (26> 

ot ^— ' ^— ' ov± vj_y — Ik ovj_ 

k (=— oo 

As usual in the weakly nonlinear theory (or weak turbulence the- 
ory), in order to close the system one must engage an equation for $1^- 
Averaging over the fast oscillations, we have 

(d/dt)\^\ 2 = 29^|cl» lik | 2 , (27) 

where suffixes k denote the kth Fourier component. 
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Equations (26) and (27) form the closed system of weakly non- 
linear equations for Jeans oscillations of the rotating homogeneous 
disk of stars, and describe a diffusion in velocity space. The spec- 
trum of oscillations and their growth rate are given by Eq. (21) and 
~ y / 2irGao\k\F(x) < Q, respectively. Thus after a time 

t ~ ~ ir 1 

the Jeans-unstable disk will assume the form described by Eq. (17). 
In the Solar vicinity, Q^ 1 ~ 2 10 8 yr. A very important feature of the 
instability under consideration is the fact that it is aperiodic (the real 
part of the wavefrequency vanishes in a rotating frame we are using) . 

As an application of the theory we investigate the relaxation of 
low frequency and Jeans-unstable, |cj*| < k and u>1 < 0, respectively, 
oscillations in the homogeneous galactic disk. Indeed, already in the 
1940s it was observed that in the Solar neighborhood the random 
velocity distribution function of stars with an age t > 10 8 yr is close to 
a Schwarzschild distribution — a set of Gaussian distributions along 
each coordinate in velocity space, i.e. close to equilibrium along each 
coordinate. In addition, older stellar populations have a higher velocity 
dispersion than younger ones. On the other hand, a simple calcula- 
tion of the relaxation time of the local disk of the Milky Way due 
to pairwise star-star encounters brings the standard value ~ 10 14 yr 
(Chandrasekhar, 1960; Binney and Tremaine, 1987, p. 489), which 
considerably exceeds the lifetime of the universe. According to our ap- 
proach, collisionless collective-type relaxation does play a determining 
role in the evolution of stellar populations of the Galactic disk. 

Evidently, the unstable fluid-like Jeans oscillations must influence 
the distribution function of the main, nonresonant part of stars in such 
a way as to hinder the wave excitation, i.e. to increase the velocity 
dispersion. This is because the Jeans instability, being essentially a 
gravitational one, tends to be stabilized by random motions (Griv and 
Peter, 1996). Therefore, along with the growth of the oscillation ampli- 
tude, random velocities must increase at the expense of circular motion, 
and finally in the disk there can be established a quasi-stationary 
distribution so that the Jeans-unstable perturbations are completely 
vanishing and only undamped Jeans-stable waves remain: Jeans insta- 
bilities are thought to heat disks up to velocity dispersion values of 

>CM-B 

2 In turn, the Jeans-stable perturbations are subject to a resonant Landau-type 
instability (Griv et at, 1999b; Griv et al, 2000b). The resonant interaction of stars 
with the field of the Jeans-stable waves will be accompanied by an additional increase 
in the velocity dispersion to values c r greater than it follows from approximate 
criterion (25). 
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In the following, we restrict ourselves to the most "dangerous," 
in the sense of the loss of gravitational stability, long-wavelength per- 
turbations, x 2 an d i 2 < 1 (see the explanation after Eq. [21]). Then 
in Eqs. (21) and (26) one can use the expansions Jf(x) ~ X 2 /4 and 
e- x h(x) « (l/2)x-(l/2)x 2 + (5/16)x 3 . Equation (26) takes the simple 
form 

f ^ 4f> < 28 > 

where D = (7r/16re 2 ) ^ k k^u^i^ 2 , and both Qa;* and $1^ are 
functions of t. As is seen, the velocity diffusion coefficient for non- 
resonant stars D is independent of v± (to lowest order). This is a 
qualitative result of the nonresonant character of the star's interaction 
with collective aggregates. 

By introducing dr/dt = D(t) and d/dt = {dr / dt){d/ dr) (Alexan- 
drov et al, 1984; Krall and Trivelpiece, 1986), Eqs. (27) and (28) are 
rewritten as follows: 

dfo d 2 fo n dD 

97-^1=°' 57 = 2 ^*' (29) 



const / v 2 



which has the solution 



f ( v± , T ) = -j= r exp\-^\. (30) 

(We have taken into account the observations that the distribution of 
newly born stars is close to the 5- function distribution, fo(v±, t = 0) = 
£( v -l); Grivnev and Fridman, 1990.) 

Two general physical conclusions can be deduced from the solution 
(30). First, during the development of the Jeans instability, Qu;* > 0, 
the Schwarzschild distribution of random velocities (a Gaussian spread 
along v r ,v<p coordinates in velocity space) is established. It sharp con- 
trast to a normal gas with frequent interparticle collisions, it is the 
collisionless collective-type interactions between the stars which trans- 
form an arbitrary initial distribution into the Maxwellian-like form. 
Secondly, the energy of the oscillation field oc ^ k l^kl 2 plays the role 
of a "temperature" T in the nonresonant-particle distribution. As the 
perturbation energy increases, the initially monoenergetic distribution 
spreads (fo becomes less peaked), and the effective temperature (or v\) 
grows with time (a Gaussian spread increases): T = 2r oc j D(t)dt oc 
/ Ylk &*3ta'*|^i,k| 2 <^- That is, during the development of the Jeans 
instability the planar mean-square velocity increases linearly with time. 
In Section 4.2 of the present paper, we verify the v \ oc t dependence 
by iV-body simulations. 
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From the above, this nonlinear wave-star interaction increases the 
velocity dispersion of stars in Milky Way's disk after they are born 
(and of particles in Jeans- unstable iV-body models). Subsequently, suf- 
ficient velocity dispersion prevents Jeans' instability from occuring. 
The "diffusion" of nonresonant stars takes place because they gain 
mechanical (oscillatory) energy as the instability develops. The ve- 
locity diffusion, however, presumably tapers off as Jeans stability is 
approached: the radial velocity dispersion c r becomes greater than the 
critical one cm- This is large enough to turn off the Jeans instability in 
a stellar disk. Thus, the true time scale for relaxation in the Milky Way 
may be much shorter than its standard value ~ 10 14 yr for the classical 
Chandrasekhar-Ogorodnikov collisional relaxation; it may be of the 
order (Qlj*) -1 > Q^ 1 > 10 9 yr, i.e. comparable with 10 periods of the 
Milky Way rotation in the Solar vicinity. The above relaxation time is 
in fair agreement with both observations and iV-body simulations of 
Milky Way's disk (Binney and Tremaine, 1987, p. 478; Grivnev and 
Fridman, 1990; Hohl, 1971). 



6. Summary 

The results show that a velocity dispersion given by Toomre's (1964) 
criterion will stabilize a disk only against axisymmetric ring-like gravity 
perturbations. However, such disks are unstable against nonaxisymmet- 
ric spiral-like perturbations. Jeans-unstable spiral perturbations "heat" 
the system; the mean-square random velocity increases linearly with 
time. The results are in agreement with analytical predictions. 

We were able to generate an axisymmetric, gravitationally stable 
disk. The initial condition for the axisymmetric stable disk was ob- 
tained analytically and confirmed experimentally: the radial dispersion 
of random velocities of stars should be equal (or greater than) to the 
modified dispersion at each radii. 
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